Sca2003-09: Multiple-point Statistics to Generate Geologically Realistic Pore-space Representations

The resolution of current micro-CT scanning –a few microns – means that for many rocks, significant porosity cannot be imaged directly in three dimensions. Furthermore, alternative approaches, such as reconstruction through simulating sedimentation and diagenesis, may be problematic for many materials whose depositional and diagenetic history is uncertain or complex. In these cases it is necessary to find another approach to generate a pore space representation. One method is to use two-point statistics to produce a three-dimensional structure from two-dimensional information obtained from thin-section analysis. However, it has been shown that these reconstructed structures often fail to reproduce the long-range connectivity of the pore space. An attractive alternative approach is to use multiple-point statistics based on twodimensional thin-sections as training images to generate geologically realistic threedimensional pore space representations that preserve the long-range connectivity of the pore structure. The method is borrowed from geostatistical techniques that use a pixelbased representation to reproduce large-scale patterns. Thin-section images can provide multiple-point statistics, which describe the statistical relation between multiple spatial locations. We test the method on a simple sandstone for which a three-dimensional image from micro-CT scanning is available. We show that the use of multiple-point statistics allows the long-range connectivity of the structure to be preserved, in contrast to two-point methods that tend to underestimate the connectivity of the pore space. The proper selection of multiple-point statistics is a key issue and is discussed in detail. INTRODUCTION Transport properties such as relative permeability and capillary pressure functions define flow behavior in porous media. These functions critically depend on the geometry and topology of the pore space, the physical relationship between rock grains and the fluids, and the conditions imposed by the flow process. A quantitative prediction of petrophysical properties in disordered media, such as sedimentary rock, frequently employs representative microscopic models of the microstructure as input. It is necessary that proper pore structural information is supplied as input to predict fluid flow properties using network models [1] or other approaches, such as the lattice-Boltzmann method [2]. The reconstruction of three-dimensional porous media is of great interest in a wide variety of fields, including petroleum engineering, biology, and medicine. An effective recons truction method allows us to generate realistic structures at will and subsequent analysis can be performed on the image in order to measure desired macroscopic transport and mechanical properties. Several methods have been proposed to generate three-dimensional pore space images. A series of two-dimensional sections can be combined to form a three-dimensional image. However, this is a delicate and laborious operation limited by the impossibility of preparing cross sections with a spacing of less than about 10μm [3]. Another approach is to use non-destructive X-ray computed microtomography [4] to image the three-dimensional pore space directly at resolutions of around a micron. However, this resolution is not sufficient to image the sub-micron size pores that are abundant in carbonates , which can only be imaged by 2D techniques such as SEM. In the absence of direct experimental 3D structural data, other reconstructions from readily available 2D microscopic images are the only viable alternative. Stochastic methods Two-dimensional thin sections are, in contrast, often readily available with high resolution. Geometrical properties such as the porosity and the two-point correlation function can be measured from these sections and used to generate a three -dimensional image with the same statistical properties. This has the advantage of being quite general, allowing a wide variety of porous media to be described [5-6]. The most widely studied stochastic reconstruction technique is based on conditioning and truncation of Gaussian random fields. This method traditionally generates a 3D binary image given its porosity and twopoint autocorrelation function. These two constraints have been found to be insufficient to reproduce the microstructure of particulate media, such as grain or sphere packs [7-9]. Including other morphological descriptors such as pore and solid phase chord distributions can improve the results [10]. As we discuss below, the disadvantage is that the simulated pore space may ignore the long-range connectivity of the pore space [11-12]. One alternative is to reconstruct the porous medium by modelling the geological process by which it was made [11-17]. New Methodology – Multiple -point Statistics While geological reconstruction overcomes the problem associated with statistical techniques, there are many systems for which a direct three-dimensional representation of the pore space is very difficult or impossible to generate. Most carbonates with intraand inter-granular porosity are good examples [18] . It would be very complex to use a processbased reconstruction method that mimics the geological history involving the sedimentation of irregular shapes, significant compaction and dissolution. In these cases it is necessary to find another approach to generate a pore space representation. One method is to use statistical techniques to produce a three-dimensional image from two-dimensional information from thin section analysis as mentioned before. Traditionally two-point statistics have been used to achieve this. However, as we have already said these images, for some cases, fail to reproduce the long-range connectivity of the pore space [11-12]. Therefore, we will pursue a multiple -point statistical technique first introduced in geostatistics to represent connected geological bodies, such as sand channels, at the field scale. Conceptually our problem is similar – we need to generate pore spaces that have a high degree of interconnectivity. In this paper we test the application of multiple-point statistics to generate a threedimensional isotropic representation of Fontainebleau sandstone. This relatively simple structure has also been analyzed by other methods and so we can test how well our method compares with other approaches. MULTIPLE-POINT STATISTICS METHOD Multiple-point statistics cannot be determined from sparse data; their inference requires a densely and regularly sampled image (called a training image) describing the geometries expected to exist in the real structure. For example, photographs of outcrops at the field scale and microscope images at the pore scale can be used as training images. The extraction of multiple-point statistics from the training image and their reproduction in a stochastic model consists of two steps: (1) borrowing mult iple-point statistics from training images and (2) pattern reproduction. References [19-20] give a detailed description of the application of multiple-point statistics to 2D field-scale reservoir description. We extend this method to apply 2D to 3D reproduction at the pore-scale. Borrowing Multiple -point Statistics from Training Images The training image of Figure 1 (thin-section image , 300×300 pixels, 7.5μm resolution/pixel) is scanned using a template t composed of nt locations ua and a central location u: t n h u u ,..., 1 = + = α α α (1) where the ha are the vectors describing the template. For our reconstruction, ha are the vectors of the square template (9×9 pixels) in Figure 1. The template is used to scan the training image and collect the pattern at each location u. The pattern (or the data event) measured by a template , for instance shown in Figure 2, is defined by { } t n h u i u i u dev ,..., 1 , ) ( ); ( ) ( = + = α α (2) Each point in the template shown in Figure 1 has a number to identify the pattern. The set of all data events scanned from the training image results in a training data set { } t j N j u dev Set ,..., 1 ), ( = = (3) where Set refers to the training data set constructed with template t. Nt is the number of different central locations of template t over the training image. The selection of the template geometry is quite important and various template sizes nt and geometries should be tried in order to reproduce better structures. Figure 1. Thin-section image and template (9×9) Figure 2. A pattern (or a data event) Pattern Reproduction In our modeling to reconstruct pore scale structures, binarized thin section images are used as training images that have only two states: pore or matrix. A detailed discussion of the training image can be found in the Results section. The multiple-point statistics are actually probabilities of occurrence of the data events. The probability of occurrence of any data event devn associated with the data template t can be inferred from the training image by counting the number c(devn) of replicates of devn found in the training image. A replicate should have the same geometric configuration and the same data values as devn. The multiple-point statistics can be identified to the proportion: ( ) ( ) n n n N dev c dev ≈ Pr where Nn is the size of the pattern. The key to this algorithm is the determination of local conditional probability distribution functions (cpdf). We need to evaluate the probability that the unknown attribute value S(u) takes any of 2 possible states – pore or matrix – given n nearest data during the reproduction at any unsampled location u. If multiple-point statistics are available, then the conditioning of S(u) to the single global data event devn can be considered, and the conditional probability can be identified to the training proportion. The cpdf is inferred directly and consistently from the training image . The multiple-point statistics, the geometrical structures in other words , are borrowed directly from the training image. This approach can theoretically apply to a three-dimensional field when three-dimensional structural information is available. Since it is difficult or impossible to measure 3D submicron scale data as mentioned above, our only alternative is to use 2D images to measure multiple-point statistics. In our case, in order to generate a three-dimensional structure from two-dimensional information, measured multiple -point sta tistics on one plane are rotated 90 degrees around each principal axis. In other words, measured statistics on the XY plane are transformed to the YZ and the XZ planes with an assumption of isotropy in orthogonal directions. In the presence of large-sca le structures, the use of a limited size template would not suffice to model the large -scale characteristics observed in the training image. The template size can be theoretically expanded to match the largest structure in the training image. However, the template size is limited by memory limitations in the numerical simulation. Therefore, an alternative approach can be introduced through a type of multigrid simulation where several different size templates with the same shapes are used to scan the training image. When large-scale structures exist in the training image, this multigrid approach is effective to measure the large -scale multiple -point statistics while requiring relatively little memory [19, 20]. To simply summarize multiple-point statistics reconstruction: a training image provides the characteristic pattern and cpdf within a designated template, then each point is reproduced using local surrounding points already reproduced and statistical information derived from a training image and a template. RESULTS In this paper multiple-point statistics have been applied to generate a pore space representation of Fontainebleau sandstone. We used four different size d templates with same 9×9 square shape , spanning 9×9, 18×18, 27×27 and 36×36 pixels. Larger scale templates can simply be expanded from the smallest-scale template shown in Figure 2, but only sample some of the pixels in the range of the template. In this multigrid system, a simulation is first performed on the largest template , which is four times larger than the smallest one – only one in 16 pixels is sampled on a 9×9 grid spanning 36×36 pixels. Once the first simulation is finished, the simulated values are assigned to the correct locations on the finer grid, and are used as conditioning data. Continuing this process, the reconstruction with multigrid system is achieved. We test how well the method reproduces the long-range connectivity of the pore space defined by the fraction of percolating cells in the structure. Experimental Sample – Fontainebleau Sandstone A three-dimensional non-destructive X-ray microtomographic image of Fontainebleau sandstone is used as the experimental data and training images for our method [21]. This sandstone serves as a basis for many rock physics experiments because of its simple structure [22]. Fontainebleau sandstone consists of monocrystalline quartz grains that have been well rounded through transport before being deposited. The sample we studied had well sorted grains of around 200μm in diameter. The sand was cemented by silica crystallizing around the grains (quartz overgrowths). Fontainebleau sandstone exhibits intergranular porosity ranging from 0.03 to roughly 0.3 depending on the degree of compaction and diagenesis. The porosity of the microtomographic image is 0.1355 compared to 0.1484 for the larger original core plug from which the imaged sample was taken. The core plug has a permeability of 1.3D and a formation factor of 22.1. The difference between the porosity of the original core and that of the microtomographic data is due to the heterogeneity of the sandstone, the difference in sample size and the resolution of microtomography. Selection of a 2D Thin-section Figure 1 shows a two-dimensional binary image of a thin section cut from a threedimensional microstructure of Fontainebleau sandstone sample, which has 300×300 pixels with 7.5μm resolution/pixel. Although the averaged porosity in the microstructure is 0.1355, the actual porosity measured from a certain 2D plane fluctuates slice by slice as shown in Figure 3. In order to reconstruct a 3D structure effectively from a 2D cut, a representative 2D image should be selected in terms of porosity. We took the slice as shown in Figure 3 by an arrow. 0.1 0.12 0.14 0.16 0.18 0.2 0 50 100 150 200 250 300 slice number po ro si ty Figure 3. Porosity fluctuation within a subregion of Fontainebleau sandstone (average porosity=0.1355, slice size=2.25×2.25 mm, slice interval=7.5μm) Conditional 3D Reconstruction Figure 4 shows a 3D multiple-point reconstruction and its 2D slice. This reconstruction took 37 hours CPU time with an Intel Xeon 1.7GHz computer. The selected 2D thin section of Fontainebleau sandstone microtomography was used as a training image to measure the multiple -point spatial relationship and a part of this 2D thin section was also used as conditioning data. The multiple-point statistics measured in a 2D plane were rotated by 90 degrees around the principal axes in order to generate a three-dimensional structure. Isotropy in the orthogonal directions is assumed. It must be noted that 3D statistical reconstruction is not constrained to reproduce a fully connected pore space. Accordingly, a small fraction of the matrix space in statistically simulated porous media may have completely surrounded by pore phase, which is unrealistic. As visual inspection shows, noise exists in the reconstructed image because this methodology uses multiple-point statistics measured from a lower dimensional image to generate a higher dimensional structure. This means there is still missing information. However, characteristic shapes of pore and matrix are reasonably reconstructed. Figure 4. Conditional 3D reconstruction of Fontainebleau sandstone (2003) and 2D cut Figure 5 shows pore-pore autocorrelation functions of the 3D microstructure. The microstructure is defined by the binary phase function (Z(r) =1 if r belongs to pore space, Z(r) =0 otherwise). The autocorrelation function RZ(u) is defined by